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Abstract 

We study temporal correlations and multifractal properties of long river discharge 
records from 41 hydrological stations around the globe. To detect long-term correla- 
tions and multifractal behaviour in the presence of trends, we apply several recently 
developed methods [detrended fluctuation analysis (DFA), wavelet analysis, and 
multifractal DFA] that can systematically detect and overcome nonstationarities in 
the data at all time scales. We find that above some crossover time that usually is 
several weeks, the daily runoffs are long-term correlated, being characterized by a 
correlation function C{s) that decays as C{s) ~ s~'^ . The exponent 7 varies from 
river to river in a wide range between 0.1 and 0.9. The power-law decay of C{s) 
corresponds to a power-law increase of the related fluctuation function F2{s) ~ s^ 
where H = 1 — 7/2. We also find that in most records, for large times, weak multi- 
fractality occurs. The Renyi exponent r(g) for q between —10 and +10 can be fitted 
to the remarkably simple form T{q) = — ln(a'^ + W) / In 2, with solely two parameters 
a and h between and 1 with a + b > 1. This type of multifractality is obtained 
from a generalization of the multiplicative cascade model. 
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1 Introduction 



The analysis of river flows has a long history. Already more than half a cen- 
tury ago Hurst found by means of his R/S analysis that annual runoff records 
from various ri vers (includin g the Nile river) exhibit "long-range statistical 



dependencies" flHurstl Il95lh . indicating that the fluctuations in water stor- 
age and runoff processes are self-similar over a wide range of time scales, 
with no single characteristic scale. Hurst's flnding is now recognized as the 
flrst e x ample for self-affine fractal behaviour in empirical time series, see e.g. 
Feder ( 19881 ). In the 1960s, the "Hurst phenomenon" was investigated on a 



broader empirical basis for ma ny other natural phenomena (jHurst et al, 
Mandelbrot and WaUij . Il96l . 



1965[ 



The scaling of the fluctuations with time is reflected by the scaling of the power 
spectrum E{f) with frequency /, E{f) ~ /"^. For stationary time series the 
exponent (3 is related to the decay of the corresponding autocorrelation func- 
tion C{s) of the runoffs (see Eq. (1)). For P between and 1, C{s) decays by a 
power law, C{s) ~ s~"', with 7 = 1-/3 being restricted to the interval between 
and 1. In this case, the mean correlation time diverges, and the system is 
regarded as long-term correlated. For /3 = 0, the runoff data are uncorrelated 
on large time scales ("white noise"). The exponents f3 and 7 can also be deter- 
mined from a fluctuation analysis, where the departures from the mean daily 
runoffs are considered as increments of a random walk process. If the runoffs 
are uncorrelated, the fluctuation function -^2(5), which is equivalent to the 
root-mean-square displacement of the random walk, increases as the square 
root of the time scale s, F2{s) ~ ^/s. For long-term correlated data, the ran- 
dom walk becomes anomalous, and F2{s) ~ s 
is related to the exponents P and 7 via (3 = 1 
data, H is identical to the classical Hurst exponent. Recently, many studies 
using these kinds of methods have dealt with scalin g properties of hydro- 
logica l records and the underlying statistics, see e. g. Loveiov and S chertzeii 
(ll99ll):lT urcotte and Greene' ('l993l ):lGuDta et al.' ('l994'):'Tes sier et all fll99f3): 
Davis et al.,n996}: Rodriguez-Tturbe and Rinalda f 1997): Pandev et al.1 (I1 998l): 



'•^ . The fluctuation exponent H 
7 = 2H — 1. For monofractal 



Matsoi|ka.s et a1 1 fennnh ilMontanari et al l (j2nnnl) :IPeters et al.1 (j2nn2h:lLivina et al 



However, the conventional methods discussed above may fail when trends are 
present in the system. Trends are systematic deviations from the average runoff 
that are caused by external processes, e.g. the construction of a water regu- 
lation device, the seasonal cycle, or a changing climate (e.g. global warming). 
Monotonous trends may lead to an overestimation of the Hurst exponent and 
thus to an underestimation of 7. It is even possible that uncorrelated data, 
under the influence of a trend, look like long-term correlated ones when us- 
ing the above analysis methods. In addition, long-term correlated data cannot 
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simply be detrended by the common technique of moving averages, since this 
method destroys the correlations on long time scales (above the window size 
used). Furthermore, it is difficult to distinguish trends from long-term corre- 
lations, because stationary long-term correlated time series exhibit persistent 
behaviour and a tendency to stay close to the momentary value. This causes 
positive or negative deviations from the average value for long periods of time 
that might look like a trend. 

In the last years, several methods such as wavelet techniques (WT) and 
detrended fluctuation analysis (DFA), have been developed that are able 
to determine long-term correlations in the presence of trends. For details 
and applications of the methods to a large number of meteorological, cli- 
matol o gical and bio logical r ecords we refer to iPeng et al. I(ll994l): iTaaau et al.l 



Bunde et al. 



fll99.5h: iBunde et al.. (2000.) : iKantelhardt et all (|2nnil ): lArneodo et al.1 ^200'^ ) 



( 2OO2I ). The methods, described in Section 2, consider fluctua- 
tions in the cumulated runoffs (often called the "profile" or "landscape" of 
the record). They differ in the way the fluctuations are determined and in the 
type of polynomial trend that is eliminated in each time window of size s. 



In this paper, we apply these detrending methods to study the scaling of the 
fluctuations F2{s) of river flows with time s. We focus on 23 runoff records 
from international river stations spread around the globe and compare the 
results with those of 18 river stations from southern Germany. We find that 
above some crossover time (typically several weeks) F2{s) scales as with 
H varying from river to river between 0.55 and 0.95 in a nonuniversal manner 
independent of the size of the basin. The lowest exponent H = 0.55 was 
obtained for rivers on permafrost ground. Our finding is not consistent with 
the hypothesis that the s caling is universal with an exponent close to 0.75 
f)Hurst et al.Lll96,4 iFedeii . 11988.) with the same power law being applicable for 
all time scales from minutes till centuries. 



The above detrending approaches, however, are not sufficient to fully char- 
acterize the complex dynamics of river flows, since they exclusively focus on 
the variance which can be regarded as the second moment -^2(5) of the full 
distribution of the fluctuations. Note that the Hurst method actually focuses 
on the ffist moment Fi{s). To further characterize a hydrological record, we 
extend the study to include all moments Fq{s). A detailed description of the 
method, which is a multifractal g eneralization of the detrended fluctuation 
analysis ( Kantelhardt et all l2002| ) an d equivalent to the W avelet Transform 



Modulus Maxima (WTMM) method (lArneodo et al. 

tion 3. Our approach differs from the multifr a.ctal approach introduced int o 



2002[ ). is given in Sec- 



hydrology by Loveioy and Schertzer (see e. g. Schertzer and Loveiov j|l987l ) ; 
Loveiov and Schertzei ( 1991 ) ; Lavallee et al. ( 19931) :IPandev et aL ( , 99^ ) that 



was based on the concept of structure functions (jFrisch and Parisi 



19851 ) and 



on the assumption of the existence of a universal cascade model. Here we per- 
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form the multifractal analysis by studying ho w the different moments of the 
fluctu ations Fq{s) scale with time s, see also iRodrieuez-Iturbe and Rinaldol 
(|l997|). We find that at large time scales, Fg{s) scales as s'^^'^\ and a sim- 
ple functional form with two parameters (a and b), h{q) = (l/q) — [ln{a'^ + 
ln(2)] describes the scaling exponent h{q) of all moments. On small 
time scales, however, a stronger multifractality is observed that may be partly 
related to the seasonal trend. The mean position of the crossover between the 
two regimes is of the order of weeks and increases with q. 



2 Correlation Analysis 



Consider a record of daily water runoff values Wi measured at a certain hy- 
drological station. The index i counts the days in the record, i = 1,2,. . . ,N. 
To eliminate the periodic seasonal trends, we concentrate on the departures 
(pi = Wi — Wi from the mean daily runoff Wi. Wi is calculated for each 
calendar date i (e.g. April 1**) by averaging over all years in the runoff se- 
ries. In addition, we checked that our actual results remained unchanged 
when also seasonal trends in the variance have been eliminated by analysing 
(j)'. = [Wi - Wi)/(wf - Wifl'^ instead of 0^. 

The runoff autocorrelation function C(s) describes, how the persistence decays 
in time. If the 0j are uncorrelated, C(s) is zero for all s. If correlations exist 
only up to a certain number of days Sx, the correlation function will vanish 
above Sy^. For long-term correlations, C(s) decays by a power-law 

C{s) = {<Pi(t),+s) ~ s-\ < 7 < 1, (1) 



where the average (. . .) is over all pairs with the same time lag s. For large 
values of s, a direct calculation of C(s) is hindered by the level of noise present 
in the finite hydrological records, and by nonstationarities in the data. There 
are several alternative methods for calculating the correlation function in the 
presence of long-term correlations, which we describe in the following sections. 



2. 1 Power Spectrum Analysis 



If the time series is stationary, we can apply standard spectral analysis tech- 
niques and calculate the power spectrum E{f) of the time series Wi as a func- 
tion of the frequency /. For long-term correlated data, we have E{f) ~ /~'^, 
where j3 is related to the correlation exponent 7 by /? = 1 — 7. This relation 
can be derived from the Wiener- Khinchin theorem. If, instead of Wi the in- 
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tegrated runoff time series z„ = Ya=i (pi is Fourier transformed, the resulting 
power spectrum scales as E{f) ~ /~^~'^- 

2.2 Standard Fluctuation Analysis (FA) 

In the standard fluctuation analysis, we consider the "runoff profile" 

n 

Zn = Y.^^^ n = l,2,...,iV, (2) 

i=l 

and study how the fluctuations of the profile, in a given time window of size 
s, increase with s. We can consider the profile as the position of a random 
walker on a linear chain after n steps. The random walker starts at the origin 
and performs, in the ith. step, a jump of length 0j to the right, if 0j is positive, 
and to the left, if 0, is negative. 

To find how the square-fluctuations of the profile scale with s, we first divide 
each record of elements into A^^ = int(A^/s) nonover lapping segments of 
size s starting from the beginning and A^^ nonoverlapping segments of size s 
starting from the end of the considered runoff series. Then we determine the 
fluctuations in each segment p. 

In the standard fluctuation analysis, we obtain the fluctuations just from the 
values of the profile at both endpoints of each segment u, F^{h',s) = [zyg — 
and average F'^iy, s) over the 2Ns subsequences to obtain the mean 
fluctuation -p2("S), 

^-^W-I^gn".^)} ■ (3) 



By definition, F2{s) can be viewed as the root-mean-square displacement of 
the random walker on the chain, after s steps. For uncorrelated (pi values, we 
obtain Fick's diffusion law F2{s) ~ s^/"^. For the relevant case of long-term 
correlations, where C{s) follows t he power-law behav iour of Eq. (1), -^2(5) 
increases by a power law (see, e.g.. iBunde et al. 



(4) 



where the fluctuation exponent H is related to the correlation exponent 7 and 
the power-spectrum exponent (3 by 

ff = 1-7/2 = (l + /3)/2. (5) 
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For power-law correlations decaying faster than 1/s, we have H = 1/2 for 
large s values, like for uncorrelated data. 



We like to note that the standard fluctuation analysis is somewhat similar 
to the resca led range analysis introduced by Hurst (for a review see, e. g., 
Feder except that it focusses on the second moment F2{s) while Hurst 



considered the first moment Fi{s). For monofractal data, H is identical to the 
Hurst exponent. 



2.3 The Detrended Fluctuation Analysis (DFA) 



There are different orders of DFA that are distinguished by the way the trends 
in the data are eliminated. In lowest order (DFAl) we determine, for each seg- 
ment z/, the best linear fit of the profile, and identify the fluctuations by the 
variance F'^{v, s) of the profile from this straight line. This way, we eliminate 
the influence of possible linear trends on scales larger than the segment. Note 
that linear trends in the profile correspond to patch- like trends in the original 
record. DFAl has been proposed originally bv lPeng et al.l (|l994^ when analyz- 
ing correlations in D NA. It can be generalized straightforwardly to eliminate 
higher order trends ( Bunde et al. . 20001 : Kantelhardt et al. . 2001 ). 



In second order DFA (DFA2) one calculates the variances F'^{h>, s) of the pro- 
file from best quadratic fits of the profile, this way eliminating the influence of 
possible linear and parabolic trends on scales larger than the segment consid- 
ered. In general, in ?7,th-order DFA, we calculate the variances of the profile 
from the best nth-order polynomial fit, this way eliminating the influence of 
possible (n — l)th-order trends on scales larger than the segment size. 

Explicitly, we calculate the best polynomial fit yu{i) of the profile in each of 
the 2Ns segments u and determine the variance 



F\u,s)^-±[z^^_ 
^ i=i 



(6) 



Then we employ Eq. (3) to determine the mean fluctuation F2{s). 

Since FA and the various stages of the DFA have different detrending capa- 
bilities, a comparison of the fluctuation functions obtained by FA and DFAn 
can yield insight into both long-term correlations and types of trends. This 
cannot be achieved by the conventional methods, like the spectral analysis. 
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2-4 Wavelet Transform (WT) 



The wavelet methods we employ here are based on the determination of the 
mean values z^{s) of the profile in each segment u (of length s), and the 
calculation of the fluctuations between neighbouring segments. The different 
order techniques we have used in analyzing runoff fluctuations differ in the 
way the fluctuations between the average profiles are treated and possible 
nonstationarities are eliminated. The first-, second- and third-order wavelet 
method are described below. 

(i) In the first-order wavelet method (WTl), one simply determines the fluctu- 
ations from the first derivative F'^^u, s) = [z,,{s) — z^+i{s)Y. WTl corresponds 
to FA where constant trends in the profile of a hydrological station are elimi- 
nated, while linear trends are not eliminated. 

(ii) In the second-order wavelet method (WT2), one determines the fluctua- 
tions from the second derivative F^(z/, s) = \z^{s) — 2z^+i{s) + ^i/+2(s)]^. So, 
if the profile consists of a trend term linear in s and a fluctuating term, the 
trend term is eliminated. Regarding trend-elimination, WT2 corresponds to 
DFAl. 

(iii) In the third-order wavelet method (WT3), one determines the fluctuations 
from the third derivative -F^(z/, s) = [zu{s) — 3ziy+i{s) + 3zt^+2{s) — ^1^+3(5)]^. 
By definition, WT3 eliminates linear and parabolic trend terms in the profile. 
In general, in WTn we determine the fluctuations from the nth derivative, 
this way eliminating trends described by (n — l)st-order polynomials in the 
data. 

Methods (i-iii) are called wavelet methods, since they can be interpreted as 
transforming the profile by discrete wavelets representing first-, second- and 
third-order cumulative derivatives of the profile. The first-order wavelets are 
known in the literature as Haar wavelets. One can also use different shapes of 
the wavelets (e . g. Ga ussian wavelets with width s), which have been used by 



Arneodo et al.l 1)2002) to study, for example, long-range correlations in DNA. 



Since the various stages of the wavelet methods WTl, WT2, WT3, etc. have 
different detrending capabilities, a comparison of their fluctuation functions 
can yield insight into both long-term correlations and types of trends. 

At the end of this section, before describing the results of the FA, DFA, and 
WT analysis, we note that for very large s values, s > N/4 for DFA and 
s > N/10 for FA and WT, the fluctuation function becomes inaccurate due 
to statistical errors. The difference in the statistics is due to the fact that the 
number of independent segments of length s is larger in DFA than in WT, 
and the fluctuations in FA are larger than in DFA. Hence, in the analysis we 
will concentrate on s values lower than Smax = N/A for DFA and Smax = N/10 
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Fig. 1. (a) The fluctuation functions F2{s) versus time scale s obtained from FA and 
DFA1-DFA5 in double logarithmic plots for daily runoff departures = ~ Wi 
from the mean daily runoff Wi for the river Weser, measured from 1823 till 1993 by 
the hydrological station Vlotho in Germany, (b) The analog curves to (a) when the 
are randomly shuffled. 

for FA and WT. When determining the scaling exponents H using Eq. (4), we 
manually chose an appropriate (shorter) fitting range of typically two orders 
of magnitude. 



2. 5 Results 



In our study we analyzed 41 runoff records, 18 of them are from the southern 
part of Germany, and the rest is from North and South America, Africa, 
Australia, Asia and Europe (see Table 1). We begin the analysis with the 
runoff record for the river Weser in the northern part of Germany, which has 
the longest record (171 years) in this study. Figure 1(a) shows the fluctuation 
functions F2{s) obtained from FA and DFA1-DFA5. In the log-log plot, the 
curves are approximately straight lines for s above 30 days, with a slope H ^ 
0.75. This result for the Weser suggests that there exists long-term persistence 
expressed by the power-law decay of the correlation function, with an exponent 
7 0.5 [see Eq. (4)]. 

To show that the slope H 0.75 is due to long-term correlations and not 
due to a broad probability distrib ution (Joseph- versus Noah- Phenomenon, 
see Mandelbrot and WallisI ( IQBSh ). we have eliminated the correlations by 



randomly shuffling the 0j. This shuffling has no effect on the probability dis- 
tribution function of 0j. Figure 1(b) shows F2{s) for the shuffled data. We 
obtain H = 1/2, showing that the exponent H ^ 0.75 is due to long-term 
correlations. 

To show that the slope if 0.75 is not an artefact of the seasonal dependence 
of the variance and skew, we also considered records where was divided by 
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10 10 
s (days) 

Fig. 2. The fluctuation functions F2{s) versus time scale s obtained from FA and 
DFA1-DFA4 in double logarithmic plots for four additional representative hydrolog- 
ical stations: (a) the Zaire in Kinshasa, Zaire, (b) the Orinoco in Puente Angostura, 
Venezuela, (c) the Mary River in Miva, Australia, and (d) the Gaula River in Haga 
Bru, Norway, (e-h) the fluctuation functions F2{s) obtained for the same rivers as 
in (a-d), from first to fifth order wavelet analysis (WT1-WT5). The straight lines 
are best linear fits to the DFAl and the WT2 results on large time scales. 

the variance of each calendary d ay and applied furth er detrending techniques 
that take into account the skew ( Livina et al. . 2003b| ). In both cases, we found 
no change in the scaling behaviour for large times (see also Sect. 3.4). This 
can be understood easily, since kinds of seasonal trends cannot effect the 
fluctuation behaviour on time scales well above one year. It is likely, however, 
that the seasonal dependencies of the variance and possibly also of the skew 
contribute to the behaviour at small times, where the slope H is much larger 
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than 0.75 in most cases (see also Sect. 3.4, where a seasonal trend in the 
variance is indeed used for modelling the crossover). 



Figure 2 shows the fluctuation functions F2{s) of 4 more rivers, from Africa, 
South America, Australia, and Europe. The panels on the left-hand side show 
the FA and DFAl-4 curves, while the panels on the right-hand side show the 
results from the analogous wavelet analysis WT1-WT5. Most curves show 
a crossover a t small time scales; a similar crossover has been reported by 
Tessier et al. for small French rivers without artificial dams or reser- 



voirs. Above the crossover time, the fluctuation functions (from DFAl-4 and 
WT2-5) show power-law behaviour, with exponents H ^ 0.95 for the Zaire, 
H ~ 0.73 for the Orinoco, H ~ 0.60 for the Mary river, and H ^ 0.55 for the 
Gaula river. 

Accordingly, there is no universal scaling behaviour since the long-term ex- 
ponents vary strongly from river to river and reflect the fact that there exist 
different mechanisms for floods where each may induce different scaling. This 
is in contrast to climate data, where universal long-term persistence of tem- 
perature records at land stations was observed (iKoscielnv-Bunde et al. 19981 : 



Talkner and Webeil . I2OOOI : IWeber and Talkneii . l200l[ lEichner et al.1 . \2Q(m 



The Mary river in Australia is rather dry in the summer and the Gaula river 
in Norway is frozen in the winter. For the Mary river, the long-term exponent 
H ~ 0.60 is well below the average value. For the Gaula river, the long-term 
correlations are not pronounced {H = 0.55) and even hard to distinguish from 
the uncorrelated case H = 0.5. We obtained similar results for the other two 
"frozen" rivers (Tana from Norway and Dvina from Russia) that we analysed. 
For interpreting this distinguished behaviour of the frozen rivers we like to 
note that on permafrost ground the lateral inflow (and hence the indirect 
contribution of the water storage in the catchment basin) contributes to the 
runoffs in a different way than on normal ground, see also Gupta and Dawdy 
[1995]. Our results (based on 3 rivers only) suggest that the contribution of 
snow melting leads to less correlated runoffs than the contribution of rainfall, 
but more comprehensive studies will be needed to confirm this interesting 
result. 

Figure 3(a) and Table 1 summarize our results for H. One can see clearly that 
the exponents H do not depe nd systematically o n the basin area A. This is 



in line with the conclusions of iGupta et al.l (|l99J) for the flood peaks, where 



a systematic dependence on A could also not be found. 

There is also no pronounced regional dependence: the rivers within a localized 
area (such as South Germany) tend to have nearly the same range of exponents 
as the international rivers. The three "frozen" rivers in our study have the 
lowest values of H. As can be seen in the figure, the exponents spread from 
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77 
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0.70 


0.53 
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0.41 


Neckar 


Horb 


65 


1 118 


0.68 


0.44 


0.75 


0.78 


Neckar 


Plochingen 


79 


3 995 


0.80 


0.49 


0.65 


0.39 


Tauber 


Bad Mergentheim 


66 


1 018 


0.80 


0.44 


0.70 


0.68 


Wertach 


Biessenhofen 


77 


450 


0.66 


0.56 


0.70 


0.31 


Wiirm 


Leutstetten 


77 


413 


0.90 


0.39 


0.66 


0.77 


Wutach 


Oberlauchringen 


85 


1 129 


0.75 


0.52 


0.67 


0.37 


Vils 


Grafenmiihle 


58 


1436 


0.61 


0.50 


0.78 


0.62 



Table 1 

Table of investigated international river basins (data from Global Runoff Data Cen- 
ter (GRDC), Koblenz, Germany) and investigated South German river basins. We 
list the river and station name, the duration of the investigated daily record, the 
size of the basin area, and the results of our analysis, H = h(2), the multifractal 
quantities a, b, and Aa. 
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Fig. 3. (a) Long-term fluctuation exponents H and (b) widths Aa of tlie f{a) spectra 
for all international records (full symbols) and all records from south Germany (open 
symbols) that we analyzed, as a function of the basin area A. Each symbol represents 
the result for one hydrological station. The dashed line in (b) is a linear fit to the 
data. 



0.55 to 0.95. Since the correlation exponent 7 is related to by 7 = 2 — 2H, 
the exponent 7 spreads from almost to almost 1, covering the whole range 
from very weak to very strong correlations. 



3 Multifractal Analysis 



3. 1 Method 



For a further characterization of hydrological records it is me aningful to extend 
Eg. (3 ) by consi dering the more g eneral fluctuation functions flBarabasi and Vicsek 
( 199ll ). see also lPavis etHI ([19961 )1 . 



2N. 



2N, 

J2 



Z(v-l)s 



(7) 



s u=l 



where the variable q can take any real value except zero. For g = 2, the 
standard fluctuation analysis is retrieved. The question is, how the fluctuation 
functions depend on q and how this dependence is related to multifractal 
features of the record. 

In general, the multifractal approach is introduced by the partition function 

Z,(s) = ^|^,,-^(,_i),r~s^(^), (8) 
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where r(g) is the Renyi scahng exponent. A record is called 'monofractal', 
when r(g) depends linearly on g; otherwise it is called multifractal. 



It is easy to verify that Zq{s) is related to Fq{s) by 
r 1 1^/^ 

F,{s) = [j;^Z,{s)j . (9) 



Accordingly, Eq. (8) implies 

Fq{s) ~ (10) 

where 

h{q) = [r{q) + l]/q. (11) 



Thus, h{q) defined in Eq. (10) is directly related to the classical multifractal 
scaling exponents T(g). 

In general, the exponent h{q) may depend on q. Since for stat ionary record s, 
h{l) is identical to the well-known Hurst exponent (see e. g. iFeder 



we will call the function h{q) the generalized Hurst exponent. For monofractal 
self-affine time series, h{q) is independent of q, since the scaling behaviour 
of the variances F^(z/, s) is identical for all segments z/, and the averaging 
procedure in Eq. (7) will give just this identical scaling behaviour for all values 
of q. If small and large fiuctuations scale differently, there will be a significant 
dependence of h{q) on q: If we consider positive values of q, the segments z/ 
with large variance F'^{v, s) (i. e. large deviations from the corresponding fit) 
will dominate the average Fq{s). Thus, for positive values of g, h{q) describes 
the scaling behaviour of the segments with large fiuctuations. Usually the 
large fiuctuations are characterized by a smaller scaling exponent h{q) for 
multifractal series. On the contrary, for negative values of g, the segments 
V with small variance F^(z/, s) will dominate the average Fq{s). Hence, for 
negative values of g, h{q) describes the scaling behaviour of the segments 
with small fiuctuations, which are usually characterized by a larger scaling 
exponent. 



In th e hydrological literature ( Rodrieuez-Iturbe and Rinaldol . lOOTtlLavallee et al 



1991 one often c onsiders the generalized mass variogram Cq{X) (see also 



Davis et al. (199 



C,(A) = (|zi+A-z,r)~A^(''). (12) 
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Comparing Eqs. (7), (10), and (12) one can verify easily that K{q) and h{q) 
are related by 



h{q)=K{q)/q. (13) 



Another way to characterize a multifractal series is the singul arity sp ectrum 



spectr 

f(a), that is related to r(q) via a Legendre transform (e.g. iFeden (|l9^ 
Rodriguez-Iturbe and Rinaldol ( 19971 )). 



a = — and f{a) = qa — T{q). (14) 
dq 



Here, a is the singularity strength or Holder exponent, while /(a) denotes the 
dimension of the subset of the series that is characterized by a. Using Eq. (11), 
we can directly relate a and /(a) to h{q), 

a = h{q) + q'E^ and f{a) = q[a - h{q)] + 1. (15) 



The strength of the multifractality of a time series can be characterized by 
the difference between the maximum and minimum values of a, Omax ~ ctmin- 
When g^^^ approaches zero for q approaching ±oo, then Aa = ttmax — «min 
is simply given by Aa = h{—oo) — h{oo). 

The multifractal analysis described above is a straightforward generalization of 
the fluctuation analysis and therefore has the same problems: (i) monotonous 
trends in the record may lead to spurious results for the fluctuation exponent 
h{q) which in turn leads to spurious results for the correlation exponent 7, 
and (ii) nonstationary behaviour characterized by exponents h{q) > 1 can- 
not be detected by the simple method since the method cannot distinguish 
between exponents > 1, and always will yield -^2(5) ~ s in this case (see 
above). To overcome these drawbacks the multifra ctal detrended fluc t uatio n 
analysi s (MF-DFA) has been introdu c ed recentlv (iK antelhardt et"aL ()2002l ). 



see also^Ko scielnv-Bunde et aL dlQQSl) : Weber and Talkner (200 D ). According 



to IKantelhardt et all (I2OO2L l2003t) . the method is as accurate as the wavelet 



methods. Thus, we have used MF-DFA for the multifractal analysis here. In 
the MF-DFA, one starts with the DFA-fluct nations F^{s) as obtained in Eq. 
(6). Then, we define in close analogy to Eqs. (3) and (7) the generalized fluc- 
tuation function, 
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Fig. 4. The multifractal fluctuation functions -Fg(s) versus time scale s obtained from 
multifractal DFA4 for two representative hydrological stations: (a) river Weser in 
Vlotho, Germany and (b) river Danube in Orsova, Romania. The curves correspond 
to different values of g = -10, -6, -4, -2, -1, -0.2, 0.2, 1, 2, 4, 6, 10 (from the top 
to the bottom) and are shifted vertically for clarity. 

Again, we can distinguish MF-DFAl, MF-DFA2, etc., according to the order 
of the polynomial fits involved. 



3.2 Multifractal Scaling Plots 



We have performed a large scale multifractal analysis on all 41 rivers. We 
found that MF-DFA2-5 yield similar results for the fluctuation function Fq{s). 
We have also cross checked the results using the W avelet Transform Modulus 
Maxima (WTMM) me thod (lArneodo et al.n2nn2h . and always find agreement 
within the error bars ( Kantelhardt et al. . boOsf T Therefore, we present here 
only the results of MF-DFA4. Figure 4(a,b) shows two representative exam- 
ples for the fluctuation functions Fq{s), for (a) the Weser river and (b) the 
Danube river. The standard fluctuation function F2{s) is plotted in full sym- 
bols. The crossover in -^2(5) that was discussed in Sect. 2.5 can be also seen 
in the other moments. The position of the crossover increases monotonously 
with decreasing q and the crossover becomes more pronounced. We are only 
interested in the asymptotic behaviour of Fg{s) at large times s. One can see 
clearly that above the crossover, the Fg{s) functions are straight lines in the 
double logarithmic plot, and the slopes increase slightly when going from high 
positive moments towards high negative moments (from the bottom to the 
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s (days) s (days) 



Fig. 5. The multifractal fluctuation functions Fq{s) obtained from multifractal DFA4 
for four additional hydrological stations: (a) Amper in Fiirstenfeldbruck, Germany, 
(b) Wertach in Biessenhofen, Germany, (c) Susquehanna in Harrisburg, USA, (d) 
Niger in Koulikoro, Mali. The g-values are identical to those used in Fig. 4. 



top). For the Weser, for example, the slope changes from 0.65 for g = 10 to 
0.9 for g = — 10 (see also Fig. 6(b)). The monotonous increase of the slopes, 
h{q), is the signature of multifractality. 



(see Figs. 
a/2 



4(c,d)), all functions 



mcrease 



When the data are shuffled 

asymptomatically as Fg(s) ~ s^/^. This indicates that the multifractality 
vanishes under shuffling. Accordingly the observed multifractality originates 
in the long-term correlations of the record a nd is not caused by singulariti es in 
the distribution of the daily runoffs (see a1so lMa,nde1brot and W^m ()l 968h ). A 
reshufflin g- resistent multif r actali ty would indicate a 'statistical' type of non- 
linearity ( Sivapalan et al. . 2002f l. We obtain similar patterns for all rivers. 
Figure 5 shows four more examples; Figs. 5(a,b) are for two rivers (Amper 
and Wertach) from southern Germany, while Figs. 5(c,d) are for Niger and 
Susquehama (Koulikoro, Mali and Harrisburg, USA). 



From the asymptotic slopes of the curves in Figs. 4(a,b) and 5(a-d), we obtain 
the generalized Hurst-exponents h{q), which are plotted in Fig. 6 (circles). 
One can see that in the whole g-range the exponents can be fitted well by the 
formula 

M,) ^ 1 - ,1T) 

q g in / 
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Fig. 6. The generalized Hurst exponents h{q) for the six representative daily runoff 
records analyzed in Figs. 4 and 5: (a) Amper in Fiirstenfeldbruck, Germany, (b) 
Weser in Vlotho, Germany, (c) Susquehanna in Harrisburg, USA, (d) Wertach in 
Biessenhofen, Germany, (e) Danube in Orsova, Romania, and (f) Niger in Koulikoro, 
Mali. The h{q) values have been determined by straight line fits of Fq{s) on large 
time scales. The error bars of the fits correspond to the size of the symbols. 



or 



X(g) = l + T(g) = l 



h[2 



(18) 



The formula can be derived from a modification of the multiplicative cascade 
model that we describe in Sect. 3.3. Here we use the formula only to show 
that the infinite number of exponents h{q) can be described by only two inde- 
pendent parameters, a and b. These two parameters can then be regarded as 
multifractal finger print for a considered river. This is particularly important 
when checking models for river flows. Again, we like to emphasize, that these 
parameters have been obtained from the asymptotic part of the generalized 
fluctuation function, and are therefore not affected by seasonal dependencies. 

We have fitted the h{q) spectra in the range — 10 < g < 10 for all 41 runoff 
series by Eq. (17). Representative examples are shown in Fig. 6. The contin- 
uous fines in Fig. 6 are obtained by best fits of h{q) (obtained from Figs. 4 
and 5 as described above) by Eq. (17). The respective parameters a and b 
are listed inside the panels of each figure. Our results for the 41 rivers are 
shown in Table 1. It is remarkable that in each single case, the q dependence 
of h{q) for positive and negative values of q can be characterized well by the 
two parameters, and all fits remain within the error bars of the h{q) values. 
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Fig. 7. The multifractal spectra f{a) for two representative runoff records (a) the 
Danube in Orsova, Romania and (b) Niger in Kouhkoro, Mali. 

From we obtain r(g) (Eq. (11)) and the singularity spectrum /(a) (Eq. (15)). 
Figure 7 shows two typical examples for the Danube and the Niger. The width 
of /(tt) taken at / = characterizes the strength of the multifractality. Since 
both widths are very different, the strength of the multifractality of river 
runoffs appears to be not universal. 

In order to characterize and to compare the strength of the multifractality for 
several time series we use as a parameter the width of the singularity spectrum 
f{a) [see Eqs. (14) and (15)] at / = 0, which corresponds to the difference 
of the maximum and the minimum value of a. In the multiplicative cascade 
model, this parameter is given by 

In a — In 6 

Aq; = . (19) 

In 2 ^ ^ 



The distribution of the Aa values we obtained from Eq. (19) is shown in 
Fig. 3(b), where we plot Aa versus the basin area. The figure shows that 
there are rivers with quite strong multifractal fluctuations, i.e. large Aa, and 
one with almost vanishing multifractality, i.e. Aa ~ 0. Two observation can 
be made from the figure: (1) There is no pronounced difference between the 
width of the distribution of the multifractality strength for the runoffs within 
the local area of southern Germany (open symbols in Fig. 3(b)) and for the 
international runoff records from all rivers around the globe (full symbols). In 
fact, without rivers Zaire and Grand River the widths would be the same. (2) 
There is a tendency towards smaller multifractality strengths at larger basin 
areas. This means that the river flows become less nonlinear with increasing 
basin area. We consider this as possible indication of river regulations that are 
more pronounced for large river basins. 

Our results for K{q) = 1 + r(g) (see Eq. (18) and Table 1) may be compared 
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with the functional form 



K{q) = (H' + l)q - ^^(g"' - q) q> 
a ' — 1 



(20) 



with the three paramete rs H', Ci, and a', that have been used by Loveioy. 



Schertzer, and coworkers (|Schertzer and ljOveiovl.ll!j»Ytl.Loveioy and 
199l[ iLavallee et all . 1199^ IXessier et ahl Il996l : IPandev et al.1 119981 ) 



oveiov ■ll987l:lLoveioy and Schertzer 



success- 
fully to describe the multifractal behaviour of rainf all and runoff records. The 
definit ion of K{q) we used in this paper is taken from Rodriguez-Iturbe and Rinaldol 
(jl997l ) and differs slightly from their definition. We like to note that Eq. (18) 
for K{q) is not only valid for positive q values, but also for negative q values. 
This feature allows us to determine numerically the full singularity spectrum 
f{a). In the analysis we focused on long time scales, excluding the crossover 
regime, and used detrending methods. We consider it as particularly inter- 
esting that only two parameters a and b or, equivalently, H and Aa, are 
sufficient to describe r(g) and K{q) for positive as well as negative q values. 
This strongly supports the idea of 'universal' multifractal behaviour of river 
runoffs as suggested (in different context) by Lovejoy and Schertzer. 

It is interesting to note that the generalized fluctuation functions we studied 
do not show any kind of multifractal phase transition at some critical value 
qo in the g-regime (—10 < q < 10) we analysed. Instead, our analysis shows 
a crossover at a specific time scale Sx (typically weeks) that weakly increases 
with decreasing moment q. In this paper, we concentrated on the large-time 
regime {s ^ Sx), where we obtained coherent multifractal behaviour and did 
not see any indication for a multifractal phase transition. But this does not 
exclude the possibility that at small scales a breakdown of mulifractality at 
a cr itical q-value may oc cur, as has been emphasized by Tessier et al. ( 1996| ) 
and lPandev et all (|l998h . 



3.3 Extended Multiplicative Cascade Model 



In the following, we like to motivate the 2-parameter formula Eq. (17) and 
show h ow it c an be obtained from the well known multifractal cascade model 
198i 



Barabasi and Vicsek . 1991; Kantelhardt et al. . 2002| ). In the 



( Feder . 

model, a record (pk of length N = 2"^™=''' is constructed recursively as fol- 
lows: In generation n = 0, the record elements are constant, i.e. 0^ = 1 for 
all = 1, . . . , A^. In the first step of the cascade (generation n = 1), the first 
half of the series is multiplied by a factor a and the second half of the series 
is multiplied by a factor b. This yields 0^ = a for = 1, . . . , N/2 and (pk = b 
for k = N/2 + 1, . . . , N. The parameters a and b are between zero and one, 
< a < 6 < 1. Note that we do not restrict the model to 6 = 1 — a as 
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is often done in the literature ( Feder . 198^. In the second step (generation 



n = 2), we apply the process of step 1 to the two subseries, yielding 0^ = 
for A; = 1, ... , N/A, 0^ = ab for k = iV/4 + 1, . . . , N/2, = ba = ab for 
k = N/2 + 1, . . . , 3iV/4, and (pk = b'^ for k = 3N/A + 1, . . . , A^. In general, in 
step n + 1, each subseries of step n is divided into two subseries of equal length, 
and the first half of the (pk is multiplied by a while the second half is multiplied 
by b. For example, in generation n = 3 the values in the eight subseries are 
a^, a^b, a?b, ab'^, a?b, ab'^ , ab'^ , b^. After nmax steps, the final generation has 
been reached, where all subseries have length 1 and no more splitting is possi- 
ble. We note that the final record can be written as (pk = a^^'''^-'^^^'^)h'^^^^^\ 
where n{k) is the number of digits 1 in the binary representation of the index 
fc, e. g. n(13) = 3, since 13 corresponds to binary 1101. 

For thi s mult i plicative cascade model, the formula for r(q) has been de rived 



i'or tni s muit i piicative cascade model, tne lormuia tor r(q) nas been ae rived 
earlier (|Fedeil . ll988l : lBarabasi and Vicsekl . EoQli iKantelhardt et al.l . l20ol . The 



result is r(g) = [— ln(a'' + b'^) + q ln(a + b)]/ In 2 or 
q g m 2 m 2 



It is easy to see that h{l) = 1 for all values of a and b. Thus, in this form 
the model is limited to cases where h{l), which is the exponent Hurst de- 
fined originally in the R/S method, is equal to one. In order to generalize 
this multifractal cascade process such that any value of h{l) is possible, we 
have subtracted the offset Ah = ln(a + 6)/ln(2) from h{q). The constant 
offset Ah corresponds to additional long-term correlations incorporated in 
the multiplicative cascade model. For generating records without this offset, 
we rescale the power spectrum. First, we fast-Fourier transform (FFT) the 
simple multiplicative cascade data into the frequency domain. Then, we mul- 
tiply all Fourier coefficients by f~^^, where / is the frequency. This way, the 
slope P of the power spectra E{f) ~ (the squares of the Fourier coeffi- 
cients) is decreased from (3 = 2h{2) — 1 = [2 ln(a + b) — ln(a^ + 6^)]/ In 2 into 
p' = 2[h{2) - Ah]-l = -ln(a2 + 62^/ln2, which is consistent with Eq. (17). 
Finally, backward FFT is employed to transform the signal back into the time 
domain. A similar Fourier filtering technique has been used by Tessier et al.l 
( 199fih when generating surrogate runoff data. 



3.4 Comparison with Model Data 



In order to see how well the extended multiplicative cascade model fits to the 
real data (for a given river), we generate the model data as follows: (i) we 
determine a and b for the given river (by best fit of Eq. (17)), (ii) we generate 
the simple multiplicative cascade model with the obtained a and b values, and 
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Fig. 8. The fluctuation functions Fq{s) obtained from tfie multifractal DFA4 for 
surrogate series generated by the extended multiplicative cascade model with pa- 
rameters a = 0.50 and b = 0.68, that correspond to the values we obtained for the 
river Wescr. The fluctuation function Fg{s) for (a) the original (f)i series and (b) 
the shuffled series are plotted versus scale s for the same values of q as in Figs. 4 
and 5. In (c) the 0, have been multiplied by 0.1 + sin^(7ri/365) before the analy- 
sis to simulate a seasonal trend. In (d) modified values of the parameters a and b 
{a = 0.26, b = 0.59) have been used on scales s < 256 to simulate the apparent 
stronger multifractality on smaller scales observed for most rivers. For the figure, 
results from 10 surrogate series of length 140 years were averaged. 

(iii) we implement the proper long-term correlations as described above. 

Figure 8(a) shows the DFA analysis of the model data with parameters a and 
b determined for the river Weser. By comparing with Fig. 4(a) we see that the 
extended model gives the correct scaling of the fluctuation functions Fq{s) on 
time scales above the crossover. By comparing Fig. 8(b) with Fig. 4(c) we see 
that the shuffled model series becomes uncorrelated without multifractahty 
similar to the shuffled data. Below the crossover, however, the model does not 
yield the observed Fq(s) in the original data. In the following we show that 
in order to obtain the proper behaviour below the crossover, either seasonal 
trends that cannot be completely eliminated from the data or a different type 
of multifractality below the crossover, represented by different values of a and 
b, have to be introduced. 



To show the effect of seasonal trends, we have multiplied the elements 0j 
of the extended cascade model by 0.1 + sin^(7ri/365), this way generating a 
seasonal trend of period 365 in the variance. Figure 8(c) shows the DFA4 
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result for the generalized fluctuation functions, which now better resembles 
the real data than Fig. 8(a). Finally, in Fig. 8(d) we show the effect of a 
different multifractality below the crossover, where different parameters a and 
h characterize this regime. The results also show better agreement with the 
real data. When comparing Figs. 8(a,c,d) with Figs. 4 and 5, it seems that the 
Danube, Amper and Wertach fit better to Fig. 8(d), i.e. suggesting different 
multifractality for short and large time scales, while the Weser, Susquehanna 
and Niger fit better to Fig. 8(c) where seasonal trends in the variance (and 
possibly in the skew) are responsible for the behaviour below the crossover. 



4 Conclusion 



In this study, we analyzed the scaling behaviour of daily runoff time series of 
18 representative rivers in southern Germany and 23 international rivers using 
both Detrended Fluctuation Analysis and wavelet techniques. In all cases we 
found that the fluctuations exhibit self-affine scaling behaviour and long-term 
persistence on time scales ranging from weeks to decades. The fluctuation 
exponent H varies from river to river in a wide range between 0.55 and 0.95, 
showing non-universal scaling behaviour. 

We also studied the multifractal properties of the runoff time series using a 
multifractal generalization of the DFA method that was crosschecked with the 
WTMM technique. We found that the multifractal spectra of all 41 records 
can be described by a 'universal' function r(g) = — ln(a'^ + 6"^)/ In 2, which 
can be obtained from a generalization of the multiplicative cascade model and 
has solely two parameters a and h or, equivalently, the fluctuation exponent 
H = i — ln(a^ + 6^)/ln4 and the width Aa = In |/ In 2 of the singularity 
spectrum. Since our function for r(g) applies also for negative q values, we 
could derive the singularity spectra /(a) from the fits. We have calculated 
and listed the values of H, a, b, and Aa for all records considered. There 
are no significant differences between their distributions for rivers in South 
Germany and for international rivers. We also found that there is no significant 
dependence of these parameters on the size of the basin area, but there is a 
slight decrease of the multifractal width Aa with increasing basin area. We 
suggest that the values of H and Aa can be regarded as 'fingerprints' for each 
station or river, which can serve as an efficient non-trivial test bed for the 
state-of-the-art precipitation-runoff models. 

Apart from the practical use of Eq. (17) with the parameters a and b that 
was derived by extending the multiplicative cascade model and that can serve 
as finger prints for the river flows, we presently are lacking a physical model 
for this behaviour. It will be interesting to see, if physically based models, 
e.g. the random tree-structure model presented in Gupta et al. ( 19961 ). can be 
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related to the multiplicative cascade model presented here. If so, this would 
give a physical explanation for how the multiplicative cascade model is able 
to simulate river flows. 



We have also investigated the origin of the multifractal scaling behaviour by 
comparison with the corresponding shuffled data. We found that the multi- 
fractality is removed by shuffling that destroys the time correlations in the 
series while the distribution of the runoff values is not altered. After shuf- 
fling, we obtain h{q) ^1/2 for all values of q, indicating monofractal be- 
haviour. Hence, our results suggest that the multifractality is not due to the 
existence of a broad, asymmetric (singular) probability density distribution 

, but due to a specific dynamical arrange- 



ment of the values in the time series, i.e. a self-similar 'clustering' of time 
patterns of values on different time scales. We believe that our results will be 
useful also to improve the understanding of extreme values (singularities) in 
the presence of multifractal long-term correlations and trends. 

Finally, for an optimal water management in a given basin, it is essential 
to know, whether an observed long-term fluctuation in discharge data is due 
to systematic variations (trends) or the results of long-term correlation. Our 
approach is also a step forward in this directions. 
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